// @(#)root/mathmore:$Id$
// Authors: L. Moneta, 12/2006
 
 /**********************************************************************
  *                                                                    *
  * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
  *                                                                    *
  * This library is free software; you can redistribute it and/or      *
  * modify it under the terms of the GNU General Public License        *
  * as published by the Free Software Foundation; either version 2     *
  * of the License, or (at your option) any later version.             *
  *                                                                    *
  * This library is distributed in the hope that it will be useful,    *
  * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
  * General Public License for more details.                           *
  *                                                                    *
  * You should have received a copy of the GNU General Public License  *
  * along with this library (see file COPYING); if not, write          *
  * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
  * 330, Boston, MA 02111-1307 USA, or contact the author.             *
  *                                                                    *
  **********************************************************************/
 
// Header file for class GSLMultiMinFunctionAdapter
//
// Generic adapter for gsl_multimin_function signature
// usable for any c++ class which defines operator( )
//
// Created by: Lorenzo Moneta  at Fri Nov 12 16:58:51 2004
//
// Last update: Fri Nov 12 16:58:51 2004


#pragma once
#include "gsl/gsl_vector.h"
#include "gsl/gsl_multifit.h"
#include <cassert>
/**
     Class for adapting any multi-dimension C++ functor class to C function pointers used by
     GSL MultiMin algorithms.
     The templated C++ function class must implement:
 
    <em> double operator( const double *  x)</em>
    and if the derivatives are required:
    <em> void Gradient( const double *   x, double * g)</em>
 
    This class defines static methods with will be used to fill the
    \a gsl_multimin_function and
    \a gsl_multimin_function_fdf structs used by GSL.
    See for examples the
    <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Providing-a-function-to-minimize.html#Providing-a-function-to-minimize">GSL online manual</A>
 
   @ingroup MultiMin
 
  */

 
  template<class UserFunc>
  struct  GSLMultiMinFunctionAdapter {
 
    static double F( const gsl_vector * x, void * p) {
 
      UserFunc * function = reinterpret_cast< UserFunc *> (p);
      // get pointer to data from gsl_vector
      return (*function)( x->data );
    }
 
 
    static void Df(  const gsl_vector * x, void * p,  gsl_vector * g) {
 
      UserFunc * function = reinterpret_cast< UserFunc *> (p);
      (*function).Gradient( x->data, g->data );
 
    }
 
    static void Fdf( const gsl_vector * x, void * p, double *f, gsl_vector * g ) {
 
      UserFunc * function = reinterpret_cast< UserFunc *> (p);
  *f  = (*function) ( x->data );
  (*function).Gradient( x-data,g->data );
 
     // (*function).FdF( x->data, f, g->data);
    }
 
  };

   typedef double ( * GSLMultiFitFPointer ) ( const gsl_vector *, void *, gsl_vector *);
   typedef void   ( * GSLMultiFitDfPointer )   ( const gsl_vector *, void *, gsl_matrix *);
   typedef void   ( * GSLMultiFitFdfPointer ) ( const gsl_vector *, void *, gsl_vector *, gsl_matrix *);
 
 
/**
   wrapper to a multi-dim function withtout  derivatives for multi-dimensional
   minimization algorithm
 
   @ingroup MultiMin
*/
 
class GSLMultiFitFunctionWrapper {
 
public:
 
   GSLMultiFitFunctionWrapper()
   {
      fFunc.f = 0;
      fFunc.df = 0;
      fFunc.fdf = 0;
      fFunc.n = 0;
      fFunc.p = 0;
      fFunc.params = 0;
#if GSL_MAJOR_VERSION > 1
      fFunc.nevalf = 0;
      fFunc.nevaldf = 0;
#endif
   }
 
 
   /// Fill gsl function structure from a C++ function iterator and size and number of residuals
   template<class FuncVector>
   void SetFunction(const FuncVector & f, unsigned int nres, unsigned int npar  ) {
      const void * p = &f;
      assert (p != 0);
      fFunc.f   = &GSLMultiFitFunctionAdapter<FuncVector >::F;
      fFunc.df  = &GSLMultiFitFunctionAdapter<FuncVector >::Df;
      fFunc.fdf = &GSLMultiFitFunctionAdapter<FuncVector >::FDf;
      fFunc.n = nres;
      fFunc.p = npar;
      fFunc.params =  const_cast<void *>(p);
   }
 
   gsl_multifit_function_fdf * GetFunc() { return &fFunc; }
 
 
  private:
 
   gsl_multifit_function_fdf fFunc;
 
};

 enum GSLMulMT{
     CFR, //gsl_multimin_fdfminimizer_conjugate_fr
     CPR, // gsl_multimin_fdfminimizer_conjugate_pr
     BFGS, //gsl_multimin_fdfminimizer_vector_bfgs 
     BFGS2, //gsl_multimin_fdfminimizer_vector_bfgs2 
     SDT // gsl_multimin_fdfminimizer_steepest_descent
 };
class GSLMultiMinimizer {
 
public:
 
   /**
      Default constructor
   */
   GSLMultiMinimizer (GSLMulMT type)  :
      fMinimizer(0),
      fType(0),
      fVec(0)
   {
      switch(type)
      {
      case CFR :
         fType = gsl_multimin_fdfminimizer_conjugate_fr;
         break;
      case CPR :
         fType = gsl_multimin_fdfminimizer_conjugate_pr;
         break;
      case BFGS :
         fType = gsl_multimin_fdfminimizer_vector_bfgs;
         break;
      case BFGS2 :
#if (GSL_MAJOR_VERSION > 1) || ((GSL_MAJOR_VERSION == 1) && (GSL_MINOR_VERSION >= 9))
         // bfgs2 is available only for v>= 1.9
         fType = gsl_multimin_fdfminimizer_vector_bfgs2;
#else
         MATH_INFO_MSG("GSLMultiMinimizer","minimizer BFSG2 does not exist with this GSL version , use BFGS");
         fType = gsl_multimin_fdfminimizer_vector_bfgs;
#endif
         break;
      case SDT:
         fType = gsl_multimin_fdfminimizer_steepest_descent;
         break;
      default:
         fType = gsl_multimin_fdfminimizer_conjugate_fr;
         break;
      }
 
   }
 
   /**
      Destructor
   */
   ~GSLMultiMinimizer ()  {
      if (fMinimizer != 0 ) gsl_multimin_fdfminimizer_free(fMinimizer);
      // can free vector (is copied inside)
      if (fVec != 0) gsl_vector_free(fVec);
   }
 
private:
   // usually copying is non trivial, so we make this unaccessible
 
   /**
      Copy constructor
   */
   GSLMultiMinimizer(const GSLMultiMinimizer &) {}
 
   /**
      Assignment operator
   */
   GSLMultiMinimizer & operator = (const GSLMultiMinimizer & rhs)  {
      if (this == &rhs) return *this;  // time saving self-test
      return *this;
   }
 
public:
 
   /**
      set the function to be minimize the initial minimizer parameters,
      step size and tolerance in the line search
    */
    template<typename FuncVector>
   int Set(const FuncVector& func, const double * x, double stepSize, double tol) {
      // create function wrapper
      fFunc.SetFunction(func);
      // create minimizer object (free previous one if already existing)
      unsigned int ndim = func.NDim();
      CreateMinimizer( ndim );
      // set initial values
      if (fVec != 0) gsl_vector_free(fVec);
      fVec = gsl_vector_alloc( ndim );
      std::copy(x,x+ndim, fVec->data);
      assert(fMinimizer != 0);
      return gsl_multimin_fdfminimizer_set(fMinimizer, fFunc.GetFunc(), fVec, stepSize, tol);
   }
 
   /// create the minimizer from the type and size
   void CreateMinimizer(unsigned int n) {
      if (fMinimizer) gsl_multimin_fdfminimizer_free(fMinimizer);
      fMinimizer = gsl_multimin_fdfminimizer_alloc(fType, n);
   }
 
   std::string Name() const {
      if (fMinimizer == 0) return "undefined";
      return std::string(gsl_multimin_fdfminimizer_name(fMinimizer) );
   }
 
   int Iterate() {
      if (fMinimizer == 0) return -1;
      return gsl_multimin_fdfminimizer_iterate(fMinimizer);
   }
 
   /// x values at the minimum
   double * X() const {
      if (fMinimizer == 0) return 0;
      gsl_vector * x =  gsl_multimin_fdfminimizer_x(fMinimizer);
      return x->data;
   }
 
   /// function value at the minimum
   double Minimum() const {
      if (fMinimizer == 0) return 0;
      return gsl_multimin_fdfminimizer_minimum(fMinimizer);
   }
 
   /// gradient value at the minimum
   double * Gradient() const {
      if (fMinimizer == 0) return 0;
      gsl_vector * g =  gsl_multimin_fdfminimizer_gradient(fMinimizer);
      return g->data;
   }
 
   /// restart minimization from current point
   int Restart() {
      if (fMinimizer == 0) return -1;
      return gsl_multimin_fdfminimizer_restart(fMinimizer);
   }
 
   /// test gradient (ask from minimizer gradient vector)
   int TestGradient(double absTol) const {
      if (fMinimizer == 0) return -1;
      gsl_vector * g =  gsl_multimin_fdfminimizer_gradient(fMinimizer);
      return gsl_multimin_test_gradient( g, absTol);
   }
 
   /// test gradient (require a vector gradient)
   int TestGradient(const double * g, double absTol) const {
      if (fVec == 0 ) return -1;
      unsigned int n = fVec->size;
      if (n == 0 ) return -1;
      std::copy(g,g+n, fVec->data);
      return gsl_multimin_test_gradient( fVec, absTol);
   }
 
 
private:
 
   gsl_multimin_fdfminimizer * fMinimizer;
   GSLMultiMinDerivFunctionWrapper fFunc;
   const gsl_multimin_fdfminimizer_type * fType;
   // cached vector to avoid re-allocating every time a new one
   gsl_vector * fVec;
 
};

//_____________________________________________________________________________________
/**
   GSLMinimizer class.
   Implementation of the ROOT::Math::Minimizer interface using the GSL multi-dimensional
   minimization algorithms.
 
   See <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html">GSL doc</A>
   from more info on the GSL minimization algorithms.
 
   The class implements the ROOT::Math::Minimizer interface and can be instantiated using the
   ROOT plugin manager (plugin name is "GSLMultiMin"). The various minimization algorithms
   (conjugatefr, conjugatepr, bfgs, etc..) can be passed as enumerations and also as a string.
   The default algorithm is conjugatefr (Fletcher-Reeves conjugate gradient algorithm).
 
   @ingroup MultiMin
*/
class GSLMinimizer  {
 
public:
 
   /**
      Default constructor
   */
   GSLMinimizer (GSLMulMT type = CFR  );
 
   /**
      Constructor with a string giving name of algorithm
    */

 
   /**
      Destructor
   */
   ~GSLMinimizer () ;
 
private:
   // usually copying is non trivial, so we make this unaccessible
 
   /**
      Copy constructor
   */
   GSLMinimizer(const GSLMinimizer &)  {}
 
   /**
      Assignment operator
   */
   GSLMinimizer & operator = (const GSLMinimizer & rhs) {
      if (this == &rhs) return *this;  // time saving self-test
      return *this;
   }
 
public:
 
   /// set the function to minimize
   template<typename FuncVector>
   void SetFunction(const  FuncVector & func);
 
 
   /// method to perform the minimization
    bool Minimize();
 
 
   /// return expected distance reached from the minimum
   double Edm() const  { return 0; } // not impl. }
 
 
   /// return pointer to gradient values at the minimum
   const double *  MinGradient() const;
 
   /// number of function calls to reach the minimum
   unsigned int NCalls() const ;
 
 
   /// minimizer provides error and error matrix
   bool ProvidesError() const { return false; }
 
   /// return errors at the minimum
   const double * Errors() const  { return nullptr; }
 
   /** return covariance matrices elements
       if the variable is fixed the matrix is zero
       The ordering of the variables is the same as in errors
   */
   double CovMatrix(unsigned int , unsigned int ) const  { return 0; }
 
 
 
 
protected:
 
private:
 
 
   GSLMultiMinimizer * fGSLMultiMin;
 
   double fLSTolerance;  // Line Search Tolerance
 
};
 
 